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Abstract —This paper proposes a fast and decentralized solution 
methodology for the robust operation of multi-area integrated 
electric-gas systems (M-IEGSs). A deterministic reformulation is 
obtained for the two-stage robust model by applying the linear- 
decision-rule based electrical reserve utilization scheme as well as 
regulating the distributed gas storages. Two linear 
approximations are developed for the nonconvex Weymouth 
equation in the gas network to determine the gas flow directions. 
The penalty convex-concave procedure (P-CCP) is then adopted to 
refine a feasible local optimum for the nonconvex model with an 
acceleration strategy. The decentralized decision-making is 
enabled by the alternating direction multipliers method (ADMM). 
The convergence as well as computational performance of the 
devised overall solution procedure can be guaranteed as only 
convex optimizations are solved. Simulation results validate the 
effectiveness of the proposed methods as well as the benefits of the 
proposed convex programing based solution procedure. 

Index Terms —Convex optimization, decentralized optimization, 
energy management, linear approximation, integrated electric-gas 
systems, robust optimization. 

Nomenclature 

Indices and Sets 

ApAf. Areas connecting to area i by tie-line/ tie-pipeline. 

c G C: Gas compressors. 

c'^/c~: Initial/terminal nodes of compressors. 

(ng): Compressors whose initial node is %. 

C~(ng): Compressors whose terminal node is ng. 

f/g G Pg: Electrical loads. 

dg^Vg. Gas loads. 

Vg{ng)\ Gas loads connected to node ng. 

Vg(s): Gas loads connected to gas storages (GSs). 

iJ^A: Areas. 

le^Cg'. Power transmission lines. 

Ip^Cp. Gas network pipelines. 

Ip Ip : Initial/terminal nodes of pipeline Ip. 

CpUg): Pipelines whose initial node is ng. 

Cp{ng)'. Pipelines whose terminal node is ng. 

^ A4’ Power network buses. 

Ug^Mg. Gas network nodes. 

s^S'. GSs. 

S{ng)'. GSs connected to node ng. 

t^T. Time periods. 


u^U^UUq: Non-gas fired units (NGUs) and Gas-fired units 
(GUs). 

uq^IAq'. GUs. 

Uq{s)\ GUs connected to GS 5*. 

Mr G zYr : Renewable generation units (RGUs). 

NGUs. 

w^W: Gas wells (GWs). 

Wing): GWs connected to node ng. 

Parameters 

Predicted values of gas loads. [MmVh] 
E7,xnsiJE^s,mm' Maximum/minimum input flow of GSs. 

[MmVh] 

^,max/^,min- Maximum/minimum output flow of GSs. 

[MmVh] 

^ax/^ni; Maximum/minimum outputs of GWs. [Mm^/h] 

Gui^: Generation shift distribution factor (GSDF) of 

NGUs and GUs. [-] 

Gu^iJGdj^'. GSDF of RGU outputs/electrical loads. [-] 

Gji^'. GSDF of tie-lines with area j. [-] 

K: Heating value of gas. [MW • h/Mm^] 

PdPPu^t'- Predicted values of electrical loads/RGU 

outputs. [MW] 

pmax/^ax; Capacities of tie-lines/tie-pipelines. [MW] 

Capacities of transmission lines. [MW] 
pmax/pmin. Capacities of NGUs and GUs. [MW] 

' Ramping rates of NGUs and GUs. [MW/h] 

Nf': Initial stored gas of GSs. [Mm^] 

^max/^min. Capacities of GSs. [Mm^] 

Zc'. Maximum compression rate of compressors. [-] 

[AFdp, AF^7]- Interval of gas loads uncertainties. [MmVh] 
[APfJ, Interval of electrical loads uncertainties. [MW] 

[AP^^, AP^py. Interval of RGUs uncertainties. [MW] 

¥i : Weymouth equation coefficients. [Mm^/(bar • 

h)2] 

^ax/^m. Nodal pressure square limits, [bar^] 

Ic' Gas consumption rates of compressors. [-] 

rjprip'. Input/output efficiency of GSs. [-] 

Gas-to-power efficiency of GUs. [-] 

Decision Variables 

fiVfsf' Input/output flow of GSs. [MmVh] 
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U' 

Flow of pipelines. [MmVh] 

fljf'a- 

Input/output flow of compressors. [MmVh] 

f : 

-'wt 

Gas production of GWs. [MmVh] 

Put- 

Outputs of NGUs and GUs. [MW] 

mV- 

Power/gas flow from area i to area j. [MW / Mm^/h] 

pij- 

Power of lines without uncertainty. [MW] 

up i^own . 

URs /DRs of NGUs and GUs. [MW] 

^st/^st • 

Stored gas in GSs. [Mm^] 

(^ut- 

RPFs of NGUs and GUs. [-] 


RPFs ofNGUs/GUs. [-] 

V’ 

Nodal pressure square, [bar^] 


Initial/terminal pressure square of pipelines, [bar^] 


Initial/terminal pressure square of compressors. 


[bar^] 


Random Variables 

Afr ' Gas loads deviations. [MmVh] 

agt 

Ap^^^ lAp^ Deviations of RGU outputs/electrical loads. [MW] 

I. Guidelines For Manuscript Preparation 

n recent years, influenced by the remarkable drop of natural 
gas prices [1], the proportion of natural gas in the total energy 

consumption rises rapidly. Compared with coal-fired units, 
natural gas-fired units (GUs) have the advantages of rapid 
dynamic characteristic and low emission, and are deployed 
gradually in power network [2]. With the growing share of GUs 
among controllable generation assets, the emerging power-to- 
gas technologies [3], and the integrated demand response [4], 
not only the physical interdependencies between power and gas 
systems have been enhanced, but also a new research topic, i.e., 
the lEGS has been emerged [5], [6]. With the broad 
participation of diverse energy suppliers, a whole integrated 
electric-gas system (lEGS) typically consists of several 
regional-IEGS (R-IEGS) owned by different operators, namely 
multi-area lEGSs (M-IEGSs), and the interconnection of 
several R-IEGS s is beneficial to the improvement of security 
and feasibility [7]. In view of the information privacy of 
operators, a decentralized optimization of M-IEGSs becomes 
an urgent need. 

Alternating direction multipliers method (ADMM) serves as 
a decentralized optimization method [8] and has been widely 
adopted in the optimization of power system [9], where the 
model convexity is required to ensure the algorithmic 
convergence. However, due to the existence of the nonconvex 
and nonlinear Weymouth equation in the gas system operation 
model, the decentralized optimization of M-IEGSs become 
challenging. Quite a few inspiring works have been carried out 
on dealing with the Weymouth equation. To provide tractable 
formulations for the commercial solvers, references [10]-[13] 
adopted the piecewise linearization method to approximate the 
Weymouth equation, transforming the original nonlinear 
program (NLP) to mixed integer linear programs (MILPs). To 
avoid additional computation burden brought by the auxiliary 
integer variables, reference [14] adopted the second-order cone 
(SOC) relaxation to tackle the Weymouth equation, where the 
gas flow directions were assumed to be fixed. Based on the 
work of [14], the penalty-convex-concave process (P-CCP) was 


adopted in [15] to drive the infeasible solution to a feasible one. 
Notably, similar solution feasibility recovery methods can be 
found in [16] and [17]. 

In practice, the gas flow directions may not be fixed, 
especially when the topology of the gas network is non-radial 
or looped. Along this line, a binary variable was introduced for 
each Weymouth equation to represent the gas flow direction, 
and the overall model would become a mixed integer second- 
order cone program (MISOCP) in [18] and [19] after 
performing the mixed integer second-order cone relaxation. To 
enhance the tightness of the relaxation, P-CCP was adopted in 
[18], where a heuristic method is proposed to improve the 
convergence of P-CCP. The big-M method is adopted in [20] 
to transform the NLP to an MISOCP, then, a multi-slack-node 
gas flow calculation method based on the Newton-Raphson 
algorithm was devised to enhance relaxation tightness. Similar 
to the SOC relaxation method, a series of linear segments were 
used to envelope SOC region through the first-order Taylor 
series expansion of Weymouth equation in references [21]-[22], 
whereas the tightness of relaxation is unconsidered. 

In addition to the internal difficulty of dealing with the 
operation model of gas systems, the uncertainties of external 
renewable power generation have also imposed great 
challenges on the operation of lEGSs, giving rise to urgent 
demands for the robust optimization method of lEGS. 
Reference [23] proposed a robust energy management 
framework for lEGSs against a set of predetermined scenarios 
of wind generation outputs. The interval-based wind generation 
uncertainty set was adopted in [24] and [25], where a robust 
dispatch strategy can be obtained for lEGSs considering 
demand response schemes. 

The aforementioned works contributed to the robust 
optimization of lEGS, whereas only a few among them focused 
on the decentralized robust optimization. The application 
scenarios of decentralized robust optimization in the existing 
works can be divided into two categories, one targeted at 
analyzing the interactions between power and gas networks 
[12], the other was designed for coordinating multiple 
interconnected R-IEGSs [18]. As mentioned above, nonconvex 
programs with binary variables (such as MILP and MISOCP) 
are used for the optimization of lEGS in the existing works, 
which challenges the convergence of ADMM. In this regard, 
reference [12] proposed a tailored ADMM, whose essence is to 
separate the process of decentralized optimization and the 
computation of binary variables. Similar to [12], reference [18] 
proposed an iterative ADMM to realize the decentralized 
optimization of M-IEGSs. In order to better understand the 
aforementioned works and provide a list of research gaps, the 
contrastive review is introduced in TABLE 1. 

TABLE I. 


THE CONTRASTIVE REVIEW OE AEOREMENTIONED WORKS 


Works 

Unfixed gas 
flow direction 

Decentralized 

Robust 

Convex 

[11][17] 

V 

X 

X 

X 

[12] 

V 

V 

V 

X 

[14][16][21] 

X 

X 

V 

V 

[15] 

X 

V 

X 

V 

[18] 

V 

V 

X 

X 

[13][20][23][24] 

V 

X 

V 

X 

[22] 

V 

X 

V 

V 
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This paper aims for the decentralized optimization of M- 
lEGSs from the perspective of eliminating nonconvexity of 
Weymouth equation, and a fast calculation method for the 
decentralized robust operation strategy of M-IEGSs is proposed. 
The interval-based uncertainty modeling as well as LDR based 
recourse action mechanism are adopted for the robust energy 
management of M-IEGSs. Two novel linear approximations for 
nonconvex Weymouth equation are derived to determine the 
directions as well as initial solutions of gas flow. Then, P-CCP 
is adopted to recover solution feasibility, where an acceleration 
method is proposed to reduce computational complexity. The 
decentralized optimization of M-IEGSs is achieved by ADMM. 
To the best knowledge of the authors, it is the first attempt of 
proposing a convex optimization based framework for the 
decentralized robust energy management for M-IEGSs. 
Compared with the existing literature, the salient features of this 
work are summarized as follows. 

(1) Two linear approximations for the nonconvex Weymouth 
equation are proposed, which are used for determining the gas 
flow directions without introducing binary variables. The 
impacts of the proposed approximation methods on the feasible 
region of the original problem are analyzed through the 
mapping between the pressure squares of the initial and 
terminal nodes of a pipeline. 

(2) A convex programming based decentralized solution 
procedure is devised for the proposed M-IEGS robust energy 
management problem, where the gas flow directions and the 
infeasibility issues in the Weymouth equations resulting from 
constraint relaxation are tackled sequentially, based on the idea 
of divide-Sind-conquer. The solutions among different regional 
lEGSs (R-IEGSs) are coordinated through the ADMM 
algorithm, whose convergence can be guaranteed by the 
convexity of the reformulated model. 

(3) Through the devised solution procedure, an optimal or 
feasible solution at least can be obtained. And the solution 
acceleration is achieved via preferentially mitigating the 
Weymouth equation infeasibility in the local gas network. 
Simulation results show that the accelerating effect could be 
more than 100 times for large-scale cases. 


II. Mathematical Eormulation 

The M-IEGSs studied includes several R-IEGSs, each R- 
lEGS is tied by power transmission lines, gas pipelines and 
communication facilities. Gas-fired units (GUs) are served as 
coupling elements for power and gas network. 


A. Objective Function 


The objective (1) is to minimize the total energy consumption 
and reserve costs of R-IEGS i. The reserve is used to mitigate 
the uncertainty of renewable generation units (RGUs) and loads. 

Cost. = rmn\y^ V [c° -fl p^.+cl ,)^] + 


Z Z K 


W _|_ ^down ^down 

U^t WxT U 




(1) 


tsT wgW 


tsT sgS 

where are the constant/linear/quadratic cost 


coefficients of NGU outputs, denotes the cost 

coefficients of upward/downward reserves of non-gas-fired 
units (NGUs), c^ is the cost coefficient of gas wells (GWs) , 
c^f express the cost coefficients of the upward/downward 
reserves of gas storages(GSs), and represents 

the flow interval of GSs. 

B. Constraints oflEGS 

In the model of lEGS, the uncertainties of RGUs and loads 
are considered as (2)-(4) [20], and the uncertainty interval can 
be determined according to the prediction accuracy. 

( 2 ) 

(3) 

AF,j<Al^,<AF,jyd^,t. (4) 

In the power network, DC power flow mode is adopted [11], 
[12], [22], and the base-case constraints without uncertainties 
are shown as (5)-(8). Particularly, (5) and (6) are the whole 
network power balancing condition and power flow constraint 
of transmission lines, respectively. Equation (7)-(8) show the 
power and ramping constraints of NGUs and gas-fired units 
(GUs). When the actual outputs of RGUs and electrical loads 
deviate from the predicted value, the outputs of NGUs and GUs 
should be adjusted accordingly. However, the power flow of 
tie-line would not be adjusted since its value has been 
determined day-ahead [26]. The operation constraints of the 
power network considering uncertainties are shown as (9)-(14). 
Particularly, the whole network power balancing condition is 
shown in (9); (10) expresses the adjustment limits of NGUs and 
GUs [27], including upward reserves (UR) and downward 
reserves (DR), whose constraints are introduced in (11). 
Equation (12)-(13) show the power and ramping constraints of 
NGUs and GUs. Equation (14) depicts the power flow 
constraint of transmission lines, where means adjustment 
of NGU and GU outputs. 

Z Pu,+ T Pu^. -Tpf-T p<i, =o> ( 5 ) 

-pr ^ pI = z G^i,Pu, + z - 


ugVI^[JV{q 

i.GAPj-i.G,jP,,<pryhd, 

jeAt 

(6) 


(7) 

1 

1 

lA 

1 

lA 

< 

(8) 

Z ^Put+^Put)-^ Z 

(9) 

-r^--<Ap„,<Cyu,t, 

(10) 

>o,c >oyu,t. 

(11) 

pr^p,,+^p^.<pryu,t. 

(12) 


(13) 
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< X (Pu. ) + Z g„_^4 ) - 

u^ 14 j^\^ 14 q u^ElIA^ 


(14) 


Z Gj,p‘i - Z G,,, {P,, +Ap,,)=p„ < Pr^h^t, 

In the gas network, the base-case operation constraints are 
shown as (15)-(23). Equations (15)-(16) describe the limits of 
gas wells (GWs) outputs and nodal pressure square. In this 
work, the steady-state gas flow model is employed, which is 
similar with [18]-[20]. Therefore, the steady-state Weymouth 
equation is presented in (17) without linepack [10], which 
depicts the relationship between the gas flow in pipeline and the 
nodal pressure squares. A simplified compressor model is 
introduced in (18)-(19) [18] [20]. The nodal gas balancing 
condition is given in (20) with linear GU model [14] [20], and 
the operation constraints of GSs are described by (21)-(23). 
Specifically, (21) is imposed to limit input and output flow of 
GSs in the worst case, (22) is added to guarantee the limit of 
available gas, and (23) sets initial stored gas for GSs. To 
mitigate the uncertainties in the gas network, the outputs of GSs 
would be adjusted. Equation (24) shows the nodal gas 
balancing condition considering uncertainties; (25)-(26) 
express the operation limits of GSs, where T means the final 
period, mean the gas flow adjustment of outputs 

and inputs of GSs. It should be noted that the uncertainty in the 
gas network mainly comes from two aspects, which are the 
inherent gas load uncertainty and the gas demand uncertainty of 
the GUs. 



^min ^ r < f 

w j Wt W ’ 9 
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> 

VI 
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> 

1 

1 
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(19) 

Z /- 
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weWirig ) 

SGSiUg) Ip^^^p(ng) 


Z j 

Z 4+ Z fc>- Z 4- 
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( 21 ) 


^s{t=T) 


IpELpiUg) ceC^iUg) ceC (tig) jsAf (Ug) 

Z Z Pu,Af?u,K = 0yn^,t, 

dg gD^ (Ug ) Uq gUq (Ug ) 

0< < F"" ,0< ,V^,L 

J St 5,max ’ d St 5,max’ ’ ’ 

r»min ^ „ , rin^in rout !^out ^ omax v-/^ > /oo\ 

(23) 

Z L+ Z (/r+A/r-r-A/i")- 

WE^iUg) SES(ng) 

Z 4-+ Z U- Z 4 

I GLp(n ) I gL (n ) cgC+( n ) 

(24) 

+ Z u- Z 4- Z 

CECirig) jeAf(ng) dgGDg(ng) 

Z (p.,,+Ap^^,)/P.,K = 0,yn^,t, 

Uq GUq {Hg ) 

<F;L,,o</r+Afr (25) 


omin ^ ~ ~ r in , \rin\^Jn 

A, ^■5^,-■5s(,-i)+(/«+A/« )'7, - 

(26) 

(fs, +A/;, )/?7. ,ys,t. 

Besides the aforementioned local operational constraints, 
capacity constraints need to be considered for the DC tie-lines 
and the gas tie-pipelines, where the gas flow in the tie-pipelines 
are assumed to be full controllable [18]. Equation (27) and (28) 
show the coordination and capacity constraints of tie-line and 
tie-pipeline, respectively. 

pf + = 0,-/p < pf < /p,V lVj g Vl (27) 
ff + ff = 0, , V/, y/ G , Vl (28) 

C. Linear Decision Rule based Uncertainty Mitigation 
Scheme 

Eor the deviations of RGUs and electrical loads with respect 
to predicted value, the output power of NGUs and GUs should 
be adjusted accordingly. Here, a Linear-decision-rule (LDR)- 
based adjustment scheme is adopted [28], as shown in (29)-(31). 
Equation (29) expresses the output power adjustments of NGUs 
and GUs are equal to Reserve participation factors (RPEs) 
multiply the total deviation of RGUs and electrical loads, RPEs 
of NGUs and GUs are continuous variables between 0 and 1 as 
equation (30). To mitigate power deviation totally, the sum of 
RPEs at each period is equal to 1 as equation (31). 


( Z Af„,, - Z Af<,,,),Vm,?. 

(29) 



0<a^^< 1, Vw,L 

(30) 

> 

II 

(31) 


ugIA'y 

Similarly, the deviations of gas loads with respect to predicted 
values need to be mitigated in real-time operation of the gas 
network, which are mitigated by the local GSs in this work [20]. 
Therefore, their input/output flow would be adjusted according 
to (32). Because GSs serve as the sole adjustable resources, the 
deviations of gas loads must be mitigated totally by local GSs, 
and the GSs must be equipped in each node connected to gas 
loads. 

AfT-Af:= Z Z I(32) 

dgEVgis) UqgUq(s) 

D. Equivalent Counterpart of the Constraints with Random 
Variables 

In constraints (5)-(28), the random variables are introduced, 
making the overall model intractable. In this regard, robust 
equivalent constraints of lEGS are proposed as (33)-(39) and 
(40)-(46) to replace (9)-(14) and (24)-(26) [29], respectively. 

Based on LDR, the power balance of whole power network 
can be satisfied. Equations (33)-(37) are imposed to preserve 
enough regulation capacity for the controllable units, which are 
minimum requirements for URs and DRs (33)-(34), base-case 
unit outputs ranges considering reserves (35), and the worst- 
case ramping constraints (36)-(37). 

rj >«„,( X AP,7 - X AP„7),V«,^. (33) 

C" ^ Z A/’™ - Z AF7), V«,r. (34) 

-C,yu,t. (35) 


- r''"""* < n < F" 

ut — rut — u 
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-CT)) ^ (36) 

(37) 


Further, (38)-(39) give the transmission line capacity 
constraints for the redistributed power flow considering the 
uncertainties as well as the active adjustments of nodal 
injections. 


pmax 


^ Pi, + Z + Z ^ (38) 


-Plt + 


°T(G„,4- Z 

,y^ 7 y I i; 


u ^ za^y lj za^y 


Guiau,)yie^Uj^,t- 


U ^ ZA'-r U ZAri, 


\u^t - ^ 

u^U^V]Uq 


..... _ 

u ^ ZAy lj ZA^y 
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u ^ zAy lj za^^ 

Je ~ ^ul^^Ut 

a^,)yi^,u^,t. 
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w ^ LJ ^^("1 


u&U^V}Ug 

- z 

u^AAy LJZ^q 


VlAet 
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>-AP^\Gjj^ 


Gul^ ^uT 

Gul^au,}yle’de,t- 


u ^ ZAy lj ZA^^ 


^ut ) 5 ’ ^* 


,>APj; 


,t y^d,i, - 2^ ^ui, 

u ^ ZAy lj 


)’ 5 ^* } 


where 6, ^,6. are instrumental variables. 

Because gas flow deviations are mitigated by the local GSs, 
the nodal gas flow balance constraint can always be satisfied. 
Robust equivalent constraints of GSs are described by (40)-(46). 
Specifically, (40)-(41) are imposed to limit input and output 
flow of GSs in the worst case, (42)-(43) are added to guarantee 
the adequacy of worst-case available gas, (44) depicts the 
temporal variation of stored gas, (45)-(46) express the 
uncertainty interval of gas loads. 

o</;<<i.wr>v^,f, (40) 

o</r (41) 

(42) 

t=\ 


Sr<s,-±Afr/rjrys,l, (43) 

t=l 

sr<s,=s,,_,,+fy:-r:‘iriT^sr,ys,t. (44) 
A/r- Z AF-- Z d-lri.K,\fs,t, (45) 

dgGVg(s) UqgUq(s) 

Afr- z z (46) 

dgGVg(s) UqgUq(s) 

where, [A/^^^,A/^^^] is the interval of aggregated gas demand 
uncertainty. 


III. Solution Methodology 

Due to the existence of the nonconvex Weymouth equation 
(17), the proposed robust energy management model for the R- 


lEGS is not readily imported to the decentralized solution 
algorithm. In what follows, the approximation and feasibility 
recovery methods for the Weymouth equation will be 
introduced first, and then the details of the adopted 
decentralized solution algorithm will be elaborated. This 
section ends with the description of the overall solution 
procedure. For the ease of analysis, the robust energy 
management problem for an R-IEGS, namely PI. 
remind) 

s.t.(5) - (8), (15) - (23), (30), (31), (33) - (46) ^ 


A. Linear Approximations for the Weymouth Equation 

According to (16)-(17), the relation between the gas flow in a 
pipeline and the pressure square difference across the pipeline, 
denoted as Aiy, is drawn in Fig. 1 and termed the Weymouth 
curve AB. In gas network, the lower and upper bounds of nodal 
pressure may be different, leading to asymmetry of the 
Weymouth curve with respect to the origin. 



Fig. 1 The geometric interpretation of the Weymouth equation and its 
approximations. 

To construct a reasonable and tractable approximation for the 
Weymouth equation, an extension of the Weymouth curve in 
the third quadrant of Fig. 1 is performed to make it a symmetric 
one as AC. Then, we connect point A and point C with a 
straight line, named the Approximation curve, and it will go 
through the origin due to the symmetry of the two points. The 
expression of the Approximation curve can be derived as 



(48) 

In this regard, we replace Weymouth equation (17) with its 
approximation (48) and preserve other constraints of PI, 
resulting in a linear approximation for PI, denoted as P2, which 
is a quadratic programming problem with linear constraints. Its 
compact form is introduced in (49), where Q , q , r , A , b are 
constant matrices and vectors, Q , q , r can be derived from 
objective function (1), x is the decision variables vectors; A ,b 
can be obtained from linear constraints (5)-(8), (15)-(16), (18)- 
(23), (30)-(31), (33)-(46) and (48). 

P2:minQ’^xQ + q^x + r 
s.t. Ax<b 

The approximation error of P2 exists inevitably, which is 
illustrated in Fig. 1. Suppose the gas flow is known, say ^ >0, 
the corresponding pressure square differences of PI and P2 can 
be obtained based on (17) and (48), and denoted as Azy and 

Azy, respectively. Intuitively, we have 
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(50) 


which can be observed from Fig. 1. In fact, a stronger 
conclusion would be drawn if//^*<0, which is 


Ar;? < Ar; 


ij 


> /;;<o, 


(51) 


considering the asymmetry of the Weymouth and the 
Approximation curves. From (50) and (51), it can be inferred 
that P2 requires larger pressure differences to maintain the same 
gas flows as PI in most occasions (two exceptions are point A 
and the origin in Fig. 1, where the pressure differences are the 
same). 

It should be noted that the pressure square difference At/^^ in 
Fig. 1 is a linear combination of two independent variables, i.e., 
and The impacts of the proposed approximated model 

on the feasible region of the nodal pressure square variables 
demonstrated in Fig. 2. 



(a) (b) 

Fig. 2 Comparison of feasible regions of PI and P2. 


Fig. 2 gathers two circumstances of the pressure mapping of 
the initial and terminal nodes with a predetermined gas flow, 
where the gas flow takes positive and negative values, 
respectively. Specially, A^+^ (^/+^) ^Fpt i^Fpt) express the 

reduced feasible intervals of the pressure square of the initial 
and terminal nodes of a pipeline, respectively, when the gas 
flow direction is positive (negative). From Fig. 2, it can be 
observed that the pressure square feasible intervals of the initial 
and terminal nodes would decrease if the Weymouth equation 
(17) is approximated by (48), regardless of the flow direction, 
leading to a smaller feasible region for the overall problem. On 
this account, if P2 has a solution, it must be a feasible (may be 
suboptimal) solution of PI. However, even if PI is feasible, P2 
can be either feasible or infeasible. 

To cope with the shrinkage of feasible region in P2, a coarser 
but still linear approximation for PI is proposed, where (16) is 
removed directly without introducing any substitute, and is 
denoted as P3, whose compact form is introduced in (52), 
where A*, b* can be obtained from linear constraints (5)-(8), 

(15) , (18)-(23), (30)-(31), (33)-(46) and (48). The only 
difference between P2 and P3 is that nodal pressure constraint 

(16) is neglected in P3. Obviously, the feasible region of P3 is 
larger than that of PI, as the values of nodal pressures are not 
constrained. 

P3: min Q^xQ + q^x + r 

. (52) 

s.t. A X <b 

From the perspective of physical meaning, the transmission 
capacity of pipeline is determined by its initial and terminal 
pressure limits. Thus, the transmission capacity of pipeline in 
PI is larger than P2, since the absolute value of pressure 


differences in P2 under the same gas flow is larger than PI [22]. 
When the limit of nodal pressure square is neglected, the 
transmission capacity of pipeline would be expanded in P3. In 
the devised overall solution procedure, P2 and P3 would be 
solved sequentially to determine the gas flow directions and 
initial gas flow, the details of which can be found at the end of 
this section. Through the combination of P2 and P3, the optimal 
or at least a feasible solution can be obtained. 


B. Feasibility Recovery 


Though P2 or P3 can provide approximate solution at 
relatively low computational costs, their solutions cannot 
reflect the actual nodal pressures of the gas network. Hence, 
additional treatment is adopted to amend the pressure 
deviations so as to recover a feasible solution of the physically 
meaningful PI. Based on the gas flow of pipeline obtained by 
P2 or P3, gas flow directions can be determined, the Weymouth 
equation can be simplified as 

(53) 


where the flow direction is assumed to be same as the positive 
direction of the pipeline for the sake of simplicity. Then, (53) 
can be replaced by a pair of opposite inequalities: 




(54) 

(55) 


Equation (54) is a convex quadratic one; according to [15], 
the concave region defined (55) can be tackled by the P-CCP 
proposed in [30] with an initial value of gas flow. The details of 
the P-CCP can be seen in Algorithm 1 [15]. Equation (56) is 
introduced to replace (55) based on the initial gas flow of 
pipeline obtained by P2 or P3, penalty term related to slack 
variable is added to objective function in (57). Thus, P4 is 
established and solved literately in Algorithm 1. 


2/;./y - (ft/ - (T,, - T ) >> 0, (56) 


where ^ f As the slack variable at iteration k. 

ipi 

P4 : min obj’f = (1) + cr^ X Z / > <^^2) 

tGTlpGCp 

S.t. (5)-(8), (15)-(16), (18)-(23), (30)-(31), (33)-(46), (54), (56) 
Note that when the flow direction is opposite, the feasibility 
recovery method is similar. 


Algorithm 1: P-CCP for feasibility recovery _ 

Step 1: Set the values of set the 

convergence tolerance S=0.l, s = 10~^, iteration index A: =0. 
Step 2: According to the solution of P2 or P3, obtain the gas 
flow directions and initial gas flows// ^ 

Step 3: Linearize the convex-concave inequality (55) with 
(56), solve P4. 

Step 4: Check the convergence residuals. 

\obj'l-obj/<8. (58) 

/ ^ //fi,, - (f/[ (59) 

If both (58) and (59) are met, then exit; else if then 

exit; else, k=k+\,Gj^=mm{KGj^_i,G^^^, set the obtained 
gas flow as the new initial valueand return to Step 3. 
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C Solution Coordination Among R-IEGSs 


The robust energy management problem for M-IEGS can be 
depicted as (60) and can be separated into sub-problems of area 
by ADMM, whose details are described in Algorithm 2 [8] [18]. 

rmn^ Costi(Xi) 
i&A 


s.t. x-GCl-,\/i, 
A^x^- +B^x^- =c^,V/,V7 


(60) 


where, x^/Costiixi) is the decision vector/ objective function of 
R-IEGS /, and defines the feasible set for x^, A^, Bg,Cg,A^, 
B^,c^ are constant matrices and vectors derived from (27)-(28). 
The augmented Lagrangian function for R-IEGS i is as follows. 


L(x^) = min 

X- gQ- 


Cojfi(x,.)+ [).J(A^Xi+B^x^.-c^) 


+^||AeX;+B^x^.-c,||^]+ Y 


2 

^yllVi + V7-‘=«ir] 




(61) 

where, and Xg are the dual variables of the coupling 
constraints, and dg are regularization parameters. 


Algorithm 2: ADMM for M-IEGSs _ 

Step 1: Initialize the value of Xy, and set convergence 
tolerance value = lO'"^, = 4x10'^, iteration index 

k=0. 

Step 2: Solve (61) for all the R-IEGSs and update and Xg 

+deiAA +B,x5 -c,). (62) 

(A^xf + B^x^ - ). (63) 

Step 3: Check the convergence residuals of each area by 


A,xf+B^xUc, 

VI 

^/.aJb^xU-x^) 

A^X; +'RgXj —Cg 

2 



where and are the tolerances of primal and dual 
residuals respectively, which can be calculated by and 
[8]. If (64) holds, terminates; otherwise, set k = k+l and 
return to Step 2._ 


D. Acceleration Strategy 

Though all the proposed models are convex ones, the potential 
computation burden could still be huge, considering the 
auxiliary constraints generated during constraint robustification, 
especially for the power flow security constraints (39). 
Therefore, an acceleration strategy is devised to ease the 
computational burden in the proposed model. 

In most cases, all the decision variables except nodal pressure 
are unchanged after amendment, which suggests the solution of 
P3 is a feasible solution of PI. Therefore, gas pressure 
deviation can be amended in gas network locally with the input 
and output of gas network fixed. The solved problem is shown 
in (65)-(73), namely P5. Based on this, the acceleration scheme 
is to replace P4 with P5 in Algorithm 1. In this case, the 


robustified power network constraints can be removed from the 
feasibility recovery process, suggesting remarkable solution 
time decrement. 

f 

P5: min of = X Z (W, + A/;,) + 

?gT weW 

Z(A/;+A/;)+ Z (A<,+A/.„;,)+Z(A//’"+A/;'-)+ 

xgS Uq gUq jeAf 

Z (Ar;;- + Ar-- + + A<r-) 

(65) 

s.t. (15)-(16), (18)-(23), (40)-(46), (54), (56), 

L=/r+V:-A/;„V>v,r, (66) 

/r - fi :=/;“■■” - /.f +a/; - a/,; , v., ?, (6?) 

= pZ + Ap.^ - Ap>,;, > Vmg , f, (68) 

fij ^ ^ g Jg ^ (69) 

C = + AC (70) 

C = C'" - AC’-. (71) 

A/:,,A/;,,A/;,A/;,Ap;,,Ap„7,AC,AC >0 (72) 

AC"’.AC;- >0 (73) 

where the symbols with prefix A and superscript III represent 
slack variables and the solution from P3, respectively. 

In some cases, the solution of P3 is infeasible with respect to 
PI and nearly all the decision variables would be changed after 
amendment. At this moment, the proposed acceleration scheme 
is ineffective, and nodal pressure deviation must be amended 
by solving P4. Now the problem is how to judge whether the 
solution of P3 is feasible for PI or not. Eirst, assuming that the 
solution of P3 is feasible for PI, then, P5 is solved by 
Algorithm 1. If the objective of P5 is beyond the threshold after 
the convergence of Algorithm 1, that means the pressure 
deviation cannot be mitigated locally in the gas network, and 
the proposed acceleration scheme is ineffective, then the P5 
would be adopted in Algorithm 1. However, if the objective of 
P5 becomes smaller than the threshold during the execution of 
Algorithm 1, that indicates the pressure deviation can be dealt 
with by the gas network and certain computation time can be 
saved. It should be noted that the proposed acceleration scheme 
would not influence the optimality of the solution as the 
variables related with objective function of original model have 
been fixed as constants. 

E. The Overall Solution Procedure 

The flowchart of the overall solution procedure is 
demonstrated in Pig. 3. Pirstly, all the R-IEGSs would 
participate in Algorithm 2 to solve P3 locally and the gas flow 
directions can be obtained. Then, the feasibility of the solution 
with respect to the unrelaxed constraints are recovered by 
solving P5 and P4 sequentially. Considering that the feasible 
region of P3 is larger than that of PI, the gas flow direction 
obtained by P3 may be infeasible, resulting the convergence 
failure of Algorithm 1&2. If so, P2 would be solved to modify 
the gas flow directions, and the feasibility of the solution would 
be recovered by solving P4. Because of the shrinkage of 
feasible region in P2, the gas flow directions obtained by P2 
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must be feasible, avoiding the convergence failure caused by 
infeasible gas flow directions. 

Noted that only convex programs are solved in Fig. 3, which 
not only provides theoretical guarantee for the convergence of 
Algorithm 1 and Algorithm 2 but also suggests acceptable 
computation burden for large-scale test systems. In order to 
clarify the differences of the five proposed models, the details 
of the five model are gathered in TABLE II, where the last two 
columns show the exclusive constraints and their descriptions. 
Concretely, PI is the original problem, which is a NLP; P2 and 
P3 are approximations of PI to obtain gas flow directions and 
initial gas flow; both P4 and P5 are used to eliminate the 
impacts of linear approximation on solution feasibility, whose 
prominent characteristic is using two convex equation (54) and 
(56) to replace the Weymouth equation. The difference of P2 
and P3 lies in the existence of nodal pressure constraints (16), 
reflecting the feasible region difference of two programs. The 
difference of P4 and P5 is whether the inputs and outputs the of 
gas network are fixed or not. 



Eig. 3 Flowchart of the overall solution proeedure 


TABLE II. 

The Comparison of the Five Proposed Programs. 



Description 

Mathematical property 

Shared constraints 

Constraints 

Description 

PI 

Original problem 

Non-linear 

programming 


(5)-(8),(33)-(39) 

(16) 

(17) 

Power network constraints 

Nodal pressure constraints 
Weymouth equation constraints 

P2 

Linear 

approximation 

Quadratic 
programming with 
linear constraints 


(5)-(8),(33)-(39) 

(16) 

(48) 

Power network constraints 

Nodal pressure constraints 

Linearized Weymouth equation constraints 

P3 

Linear 

approximation 

Quadratic 
programming with 
linear constraints 

(15), (18)-(23), 

(5)-(8),(33)-(39) 

(48) 

Power network constraints 

Linearized Weymouth equation constraints 




(30)-(31), 

(5)-(8),(33)-(39) 

Power network constraints 


Feasibility 

Quadratic 


(16) 

(54) 

(56) 

Nodal pressure constraints 

SOC constraints of Weymouth equation 
Constraints added by P-CCP 

P4 

recovery of linear 

programming with 


(16) 

Nodal pressure constraints 


approximation 

SOC constraints 


(54) 

(56) 

(66)-(73) 

SOC constraints of Weymouth equation 
Constraints added by P-CCP 
Constraints used to fix the inputs and 
outputs of gas network 



Fig. 4 Topology of the two-area M-IEGS. 


In this section, the performances of the proposed linear 
approximations as well as the acceleration scheme are justified 
by several different gas networks and single R-IEGS of 
different scales; then, the proposed methods are first 
implemented on a small-scale M-IEGS to validate the 
effectiveness, whose topology is depicted in Fig. 4; at last, the 
results of scalability tests for the proposed methods are 
demonstrated. All the parameters and topologies of all test 
systems in detail can be seen in [31]. 

The test system depicted in Fig. 4 consists of two R-IEGSs, 
denoted as R-IEGS I and R-IEGS II, respectively. R-IEGS I 
includes a 5-bus power network and a 7-node gas network, and 
R-IEGS II incorporates a 6-bus power network and a 6-node 
gas network. The two R-IEGSs are interconnected through a 
DC tie-line and a gas tie-pipeline. 

A. Performances of the LP-based Approximation Methods 

To test the performances of the proposed LP-based 
approximation methods for the Weymouth equations, the 











































































AUTHOR et al: TITLE IS LIMITED TO 50 WORDS 


9 


results of P2 and P3 under three gas networks are compared 
with the results from PI, which are solved by the IPOPT solver 
of OPTI toolbox. The results are demonstrated in TABLE III, 
where G7, GIO, and G20 are 7-, 10-, and 20-node gas networks, 
respectively. The topologies as well as related parameters of the 
test systems can be found in [31]. Specifically, the simulated 
models only incorporate operational constraints of the gas 
network so as to highlight the impacts of approximations on the 
Weymouth equations. 

Noted that the solutions provided by IPOPT are always 
feasible in the three tests. In this regard, the last column of 
TABLE III is regarded as the benchmark. Erom TABLE III, it 
can be seen that the performances of the three programs are the 
same for G7 and GIO, in terms of optimality and feasibility. In 
the G20 case, P2 offers a 2.14% larger objective than PI, while 
the solution feasibility with respect to the Weymouth equation 
is reached, which confirms P2 would cut off part of the feasible 
region of PI and cause the increase of the objective. It should 
be observed that the solution quality of P3 is the same with PI 
in the G20 case, which indicates P3 is a necessary supplement 
to P2. Noted that the P3 is solved ahead of P2 in Algorithm I, 
therefore, the proposed models can offer feasible solutions with 
good quality in all these three cases. 

TABLE III. 

Simulation Results of P1-P3 


P3 P2 PI 


Test 

system 

Objective 

function 

($) 

Sum of 
slack 
variables 

Objective 

function 

($) 

Sum of 
slack 
variables 

Objective 

function 

($) 

G7 

6.2632x10“ 

7.2943x10''“ 6.2632x10“ 

7.2943x10''“ 

6.2632x10' 

GIO 

1.4569x10'^ 

3.4835x10''’ 

1.4569x10" 

3.6473x10'" 

1.4569x10' 

G20 

1.3949x10’ 

5.7857x10'’ 

1.4247x10’ 

1.9079x10'’ 

1.3949x10' 


B. Effectiveness of the Acceleration Technique 
To validate the effectiveness of the devised acceleration 
technique in Section IV.B, i.e. solving P5 instead of P4 using 
Algorithm 1 in the middle of the overall solution procedure, the 
computation performance of Algorithm 1 with P5 and 
Algorithm 1 with P4 under 4 different R-IEGSs are gathered in 
TABLE IV. The first column of TABLE IV indicates the size 
of the test system, such as P5-G7 representing a 5-bus power 
network integrated with a 7-node gas network. In each test 
system, the intervals of the RGU outputs/electrical loads/gas 
loads uncertainty are 1/1/0.1% of the predicted values, 
respectively. The topologies as well as related parameters of the 
test systems can be found in [31]. 


TABLE IV. 

Computation Performance of Algorithm I Under P4 and P5 


Test 

System 


P5 

P4 

Number of 
reduced 
constraints 

Time 

(s) 

Iteration 

Time 

(s) 

Iteration 

P5-G7 

0.245 

1 

0.334 

1 

3648 

P6-G6 

0.413 

3 

0.997 

3 

4152 

P24-G10 

3.791 

1 

1632.9 

1 

83712 

P39-G20 

1.972 

1 

365.2 

2 

100224 


C. Strategy Robustness Test 


In the proposed model, all the robust constraints are 
transformed into normal ones. To validate the effectiveness of 
the transformed model, one thousand scenarios are generated 


randomly. The intervals of the RUs/electrical loads power/ gas 
loads uncertainty are 15/10/1.5% of the predicted values. 


4()U 200 
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(a) 





(b) 




Fig. 5 The outputs of units and power flows in random eases. 


In power network, only one GU {uQiin R-IEGS I) and one 
line (/g 3 in R-IEGS I) are tested, and the result is shown in Eig. 
5. The and indicate the power of GU in random cases 

and the base case, p. . and p. . indicate the power of line in 

IqI IqI 

random cases and the base case, respectively. Erom Eig. 5, it 
can be observed that the power adjustment of GU fall between 
UR and DR, and the power flow is always below the capacity 
of the transmission line. 

In gas network, only one GS (equipped in dgxof R-IEGS I) 
is tested, and the result is shown in Eig. 6. The s^t and s^^ 
indicate the reserves of GS in random cases and base case, ? 
and^^ indicate the flow of GS in random cases and the base 
case, where positive value means gas flow out GS and negative 
value means gas flow into GS, A/^^^and indicate the 
maximum and minimum flow of GS under worst cases, S^^^ 
and S“^^ indicate the maximum and minimum reserves of GS, 
respectively. Similarly, it can be viewed that the flow and 
adjustment of GS fall between the maximum and minimum in 
all the cases. 


U.U2 0.4 



-0.04 ---------- 0 —' 

0 5 10 15 20 25 0 5 10 15 20 25 

Timc/h Time/h 

(a) (b) 

Fig. 6 Flow and reserves of GS in random eases. 

D. Comparison Among the Different Gas Network Modeling 
Methods 

To highlight the computation benefits of convex 
programming, some other solution methods for the 
decentralized robust energy management of M-IEGSs are 
implemented, including the IPOPT method (solving the 
nonlinear and nonconvex PI directly using the OPTI toolbox). 
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methods in [18] and [10]. The simulation results are 
summarized in TABLE V, where the term “Gap” indicates 
difference of the left and right hand side terms of the Weymouth 
equation (17), and MILPs is short for mixed integer linear 
programs. It should be noted that the solution coordination 
among R-IEGSs is achieved by ADMM in all the listed 
methods, and the differences lie in the solution methods for the 
gas network models. Meanwhile, the intervals of the RGU 
outputs/electrical loads/gas loads uncertainty are 1/1/1%. 

Erom TABLE V, it can be observed that the proposed method 
can provide a good solution from the optimality and feasibility 
perspectives, in the shortest time. Among those, the 
performance of IPOPT is the worst, as it fails to generate a 
solution when r=12. Moreover, the solutions of [10] cannot 
meet the feasibility requirement in three tests, as a trade-off 
between computation burden and solution accuracy. When 
r=12, the proposed methods is almost 6 times faster than [18] 
and the solution quality is even better in terms of feasibility. 
Another observation is that the increment of computation time 
of the proposed methods is the lowest among all these methods. 
The reason is that the computation burden of the method in [18] 
and [10] grows fast along with the growth of the periods number, 
as the number of binary variables would grow accordingly, 
indicating the worse performance of the IPOPT and the 
methods in [18] and [10] in larger test systems as well as the 
necessity and significance of developing convex optimization 
based solution procedures. 

TABLE V. 


Comparison of Gas Network Modeling Methods 


Period 

(h) 

Method 

Model 

property 

Objective 

($) 

Gap (%) 

Time 

(s) 


IPOPT 

NLP 

1.7090x10*^ 

9.9920x10''^ 

1.681 

9 

Ref. [18] 

MISOCPs 

1.7090x10*^ 

3.2758x10''^ 

6.407 


Ref. [10] 

MILPs 

1.7090x10*^ 

1.48 

2.759 


Proposed 

SOCPs 

1.7090x10'^ 

4.3874x10'^ 

1.231 


IPOPT 

NLP 

5.5485x10'^ 

9.9976x10’'^ 

2.171 

6 

Ref. [18] 

MISOCPs 

5.5485x10'^ 

7.5831x10''^ 

6.810 

Ref. [10] 

MILPs 

5.5485x10*^ 

1.48 

54.823 


Proposed 

SOCPs 

5.5485x10*^ 

3.9869x10'^ 

1.822 


IPOPT 

NLP 

- 

- 

- 

12 

Ref. [18] 

MISOCPs 

1.0300x10’ 

1.2252x10'’ 

12.066 

Ref. [10] 

MILPs 

1.0300x10’ 

1.48 

217.058 


Proposed 

SOCPs 

1.0300x10’ 

1.7436x10'^ 

1.845 


E. Computational Efficiency Test 

In the sequel, the proposed methods are performed on two 
larger test systems, denoted as TS-I and TS-II, whose 
parameters can be found in [31], so as to demonstrate the 
computational performance of proposed method on large scale 
test system. The components of TS-I and TS-II are listed in 
TABLE VI, where E118G48 means the R-IEGS is composed 
of a 118-bus power network and a 48-node gas network, and so 
on in a similar fashion. Here, times of computation burden 
(TCB) is defined as the ratio of the total number of constraints 
of the test system to that of the test system in Section IV. A with 
24 time periods, which can reflect the computation burden 
increment. The simulation results are gathered in TABLE VI, 
where the column index ATPI is short for average time per 
iteration. 


Erom TABLE VI, it can be observed that the total 
computation time as well as ATPI grows almost linearly along 
with the number of periods in both test systems. Eor the day- 
ahead dispatch of the two test systems, the simulation time is 
within one hour and does not go up remarkably along with 
system scale, which reveals the benefits of the proposed 
convex-programming-based solution procedure. 

TABLE VI. 


Computational Performance on Large-Scale Test Systems 



Systems 

Time 

TCB 

Iteration 

Total 

ATPI 


periods 

time (s) 

(S) 

Bench 

E5G7- 

24 

1 


1.859 

0.266 

7 

mark 

E6G6 




4 

14 

5 

137.5 

27.5 


E39G20- 

E118G48 

8 

29 

5 

263.1 

52.6 

TS-I 

12 

57 

5 

410.1 

82.0 


16 

86 

5 

559.4 

111.9 



24 

115 

5 

1012.4 

202.5 



4 

27 

12 

292.4 

24.4 


E24G10- 

8 

55 

12 

517.9 

43.2 

TS-II 

E39G20- 

12 

82 

13 

1065.1 

81.9 


E118G48 

16 

110 

16 

2390.5 

183.9 



24 

165 

15 

3529.6 

235.3 


V. Conclusion 

In this paper, a decentralized solution procedure for the robust 
energy management of M-IEGSs is proposed to mitigate the 
uncertainties from RGU outputs and energy demands. With the 
help of the LDR-based electrical reserve utilization scheme as 
well the distributed GSs, a single-level deterministic 
equivalence can be obtained for the non-trivial robust problem. 
Novel linear approximation method is proposed to tackle the 
nonconvexities brought by the Weymouth equations in gas 
networks, and the P-CCP are employed to recover feasibility of 
approximate solution with an acceleration strategy. The 
ADMM is adopted to realize a decentralized computation of the 
robust energy management strategy. The overall solution 
procedure only requires solving convex optimizations, which 
provides convergence guarantee for the ADMM and P-CCP 
algorithms, and indicates fair computational performances, 
especially for moderate test systems. The effectiveness of the 
linear approximations is verified by simulation results. It also 
can be observed that the proposed methods can provide an 
acceptable solution under the given time limit especially for 
large-scale cases. The proposed models and methods are also 
suitable for addressing the locational marginal price based 
market clearing problems of the M-IEGS, where convex 
programs are preferred. 
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